require(stringi)
require(coda)
require(RColorBrewer)
color <- brewer.pal(9, "Set1")
paired <- brewer.pal(12, "Paired")

# load post-processed draws and produce figures and summary statistics for the paper
# theta: standardized party ideal points; beta: rescaled committee loadings
load("Wordshoal/postprocessing_out.Rdata")

# load bootstrap array to recover party ordering and recompute observation periods
load("Wordfish_bootstrap/bootstrap_array.Rdata")

n <- max(party.position.data$party.id)  # number of parties
m <- sum(stri_detect_regex(colnames(party.position.data), "[一-龠]+"))  # number of committees

party.begin <- party.end <- rep(NA, n)  # first and last observed year for each party
for (i in 1:n) {
  party.begin[i] <- min(party.position.data$time[party.position.data$party.id == i])
  party.end[i] <- max(party.position.data$time[party.position.data$party.id == i])
}
committee.begin <- committee.end <- rep(NA, m)  # first and last observed year for each committee
for (j in 1:m) {
  temp <- subset(party.position.data, ! is.na(party.position.data[, j + 3]))
  committee.begin[j] <- min(temp$time)
  committee.end[j] <- max(temp$time)
}

#### party positions ####
posterior.draws.party <- postprocessing.out[, stri_detect_regex(colnames(postprocessing.out), "theta")]

# posterior means and 95% HPD intervals for party-year ideal points
point.estimates.party <- colMeans(posterior.draws.party)
credible.intervals.party <- matrix(NA, ncol(posterior.draws.party), 2)
rownames(credible.intervals.party) <- colnames(posterior.draws.party)
for (i in 1:ncol(posterior.draws.party)) {
  credible.intervals.party[i, ] <- HPDinterval(mcmc(posterior.draws.party[, i]))
}

party.levels <- levels(factor(party.position.data$party))
party.names <- unique(stri_replace_all_regex(party.levels, "\\_.+", ""))

# Note: Party IDs correspond to the factor levels in party.position.data$party.
party.labels <- party.names
party.labels[5] <- "DPJ/DP"
party.labels[10] <- "JRP/JIP"
party.labels[11] <- "JSP/SDP"
party.labels[14] <- "Mushozoku no Kai"
party.labels[20] <- "POH"
party.labels[26] <- "SDF"

LDP_HOR.id <- 22
JSP_HOR.id <- 18
DPJ_HOR.id <- 6

major.party.HOR.id <- c(6, 10, 12, 18, 20, 22)
major.party.HOC.id <- major.party.HOR.id - 1
color.id <- c(3, 7, 4, 2, 5, 1)
lty.id <- c(6, 4, 3, 2, 5, 1)

# Figure 2: ideal-point trajectories for major parties in each chamber (means and 95% HPD intervals)
png("Figure_2.png", 
    width = 5, height = 8, unit = "in", pointsize = 12, res = 2000)
layout(matrix(1:3, 3, 1, byrow = TRUE), heights = c(1, 1, 0.25))
par(mar = c(4, 4, 3, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(1, 61), ylim = c(-6, 2), 
     main = "House of Representatives", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(-6, 2, 2), lty = 3, col = "lightgray")
for (i in 1:6) {
  polygon(c(party.begin[major.party.HOR.id[i]]:party.end[major.party.HOR.id[i]], 
            party.end[major.party.HOR.id[i]]:party.begin[major.party.HOR.id[i]]), 
          c(credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", major.party.HOR.id[i], "\\,")), 1], 
            rev(credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", major.party.HOR.id[i], "\\,")), 2])), 
          border = NA, col = paste0(color[color.id[i]], "30"))
  lines(party.begin[major.party.HOR.id[i]]:party.end[major.party.HOR.id[i]], 
        point.estimates.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", major.party.HOR.id[i], "\\,"))], 
        lty = lty.id[i], lwd = 1, col = color[color.id[i]])
}
axis(1, at = seq(1, 61, 10), labels = NA, lwd = 0.5)
mtext(c(seq(59, 99, 10), "09", "19"), 1, 1, at = seq(1, 61, 10), cex = 0.8)
mtext("Year", 1, 2.5, cex = 0.8)
axis(2, at = seq(-6, 2, 2), labels = NA, lwd = 0.5)
mtext(c("-6", "-4", "-2", "0", "4"), 2, 1, at = seq(-6, 2, 2), cex = 0.8)
mtext("Ideal point", 2, 2.5, cex = 0.8)
plot(NULL, NULL, type = "n", xlim = c(1, 61), ylim = c(-6, 2), 
     main = "House of Councillors", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(-6, 2, 2), lty = 3, col = "lightgray")
for (i in 1:6) {
  polygon(c(party.begin[major.party.HOC.id[i]]:party.end[major.party.HOC.id[i]], 
            party.end[major.party.HOC.id[i]]:party.begin[major.party.HOC.id[i]]), 
          c(credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", major.party.HOC.id[i], "\\,")), 1], 
            rev(credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", major.party.HOC.id[i], "\\,")), 2])), 
          border = NA, col = paste0(color[color.id[i]], "30"))
  lines(party.begin[major.party.HOC.id[i]]:party.end[major.party.HOC.id[i]], 
        point.estimates.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", major.party.HOC.id[i], "\\,"))], 
        lty = lty.id[i], lwd = 1, col = color[color.id[i]])
}
axis(1, at = seq(1, 61, 10), labels = NA, lwd = 0.5)
mtext(c(seq(59, 99, 10), "09", "19"), 1, 1, at = seq(1, 61, 10), cex = 0.8)
mtext("Year", 1, 2.5, cex = 0.8)
axis(2, at = seq(-6, 2, 2), labels = NA, lwd = 0.5)
mtext(c("-6", "-4", "-2", "0", "4"), 2, 1, at = seq(-6, 2, 2), cex = 0.8)
mtext("Ideal point", 2, 2.5, cex = 0.8)
par(mar = c(0, 0, 0, 0))
plot.new()
legend("center", 
       legend = c("Liberal Democratic Party", "Japanese Communist Party", 
                  "Komeito", "Japanese Socialist Party/\nSocial Democratic Party", 
                  "Democratic Socialist Party", 
                  "Democratic Party of Japan/\nDemocratic Party"), 
       lty = c(1, 3, 5, 2, 4, 6), lwd = 1, 
       col = color[c(1, 4, 5, 2, 8, 3)], bty = "n", ncol = 2)
dev.off()

# Figure A3: party-specific trajectories (means and 95% HPD intervals), with reference lines for key parties
png("Figure_A3.png", 
    width = 6, height = 8, unit = "in", pointsize = 8, res = 2000)
layout(matrix(1:28, 7, 4, byrow = TRUE))
par(mar = c(1.5, 1.5, 2, 0.2), lwd = 0.5)
for (i in 1:length(party.names)) {
  plot(NULL, NULL, type = "n", xlim = c(1, 61), ylim = c(-5, 4), 
       main = party.labels[i], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  HOR.id <- which(party.levels == paste0(party.names[i], "_HOR"))
  HOC.id <- which(party.levels == paste0(party.names[i], "_HOC"))
  if (length(HOC.id) == 1) {
    polygon(c(party.begin[HOC.id]:party.end[HOC.id], party.end[HOC.id]:party.begin[HOC.id]), 
            c(credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOC.id, "\\,")), 1], 
              rev(credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOC.id, "\\,")), 2])), 
            border = NA, col = paste0(paired[5], "30"))
  }
  if (length(HOR.id) == 1) {
    polygon(c(party.begin[HOR.id]:party.end[HOR.id], party.end[HOR.id]:party.begin[HOR.id]), 
            c(credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOR.id, "\\,")), 1], 
              rev(credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOR.id, "\\,")), 2])), 
            border = NA, col = paste0(paired[1], "30"))
  }
  lines(party.begin[LDP_HOR.id]:party.end[LDP_HOR.id], 
        point.estimates.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", LDP_HOR.id, "\\,"))], 
        col = "gray")
  lines(party.begin[JSP_HOR.id]:party.end[JSP_HOR.id], 
        point.estimates.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", JSP_HOR.id, "\\,"))], 
        col = "gray")
  lines(party.begin[DPJ_HOR.id]:party.end[DPJ_HOR.id], 
        point.estimates.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", DPJ_HOR.id, "\\,"))], 
        col = "gray")
  if (length(HOC.id) == 1) {
    lines(party.begin[HOC.id]:party.end[HOC.id], 
          point.estimates.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOC.id, "\\,"))], 
          lty = 5, col = paired[6])
  }
  if (length(HOR.id) == 1) {
    lines(party.begin[HOR.id]:party.end[HOR.id], 
          point.estimates.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOR.id, "\\,"))], 
          col = paired[2])
  }
  axis(1, at = seq(1, 61, 10), labels = NA, lwd = 0.5)
  mtext(c(seq(59, 99, 10), "09", "19"), 1, 0.5, at = seq(1, 61, 10), cex = 0.6)
  axis(2, at = seq(-5, 4, 1), labels = NA, lwd = 0.5)
  mtext(c("-4", "0", "4"), 2, 0.5, at = c(-4, 0, 4), cex = 0.6)
}
dev.off()

# reshape party-year draws into a named list for convenience
# each element is a draws-by-year matrix for one party (e.g., "LDP_HOR")
posterior.draws.party.list <- list()
for (i in 1:length(party.levels)) {
  posterior.draws.party.list[[i]] <- posterior.draws.party[, stri_detect_regex(colnames(posterior.draws.party), paste0("theta\\[", i, "\\,"))]
  rownames(posterior.draws.party.list[[i]]) <- NULL
  colnames(posterior.draws.party.list[[i]]) <- seq(party.begin[i], party.end[i]) + 1958
  names(posterior.draws.party.list)[[i]] <- party.levels[i]
}

# probability that the LDP is to the left of the JRP/JIP in 2019 (HOR)
round(mean(posterior.draws.party.list[["LDP_HOR"]][, "2019"] < posterior.draws.party.list[["JRP_HOR"]][, "2019"]), 3)

# probability that the LDP is to the left of the PJK in 2016 (HOC)
round(mean(posterior.draws.party.list[["LDP_HOC"]][, "2016"] < posterior.draws.party.list[["PJK_HOC"]][, "2016"]), 3)

# probability that the DPJ is to the left of the LDP in 2012 (HOR)
round(mean(posterior.draws.party.list[["DPJ_HOR"]][, "2012"] < posterior.draws.party.list[["LDP_HOR"]][, "2012"]), 3)

# probability that the DPJ is to the left of the JCP in 2012 (HOR)
round(mean(posterior.draws.party.list[["DPJ_HOR"]][, "2012"] < posterior.draws.party.list[["JCP_HOR"]][, "2012"]), 3)

# average distance between the DSP and the JSP (HOR)
round(mean(posterior.draws.party.list[["DSP_HOR"]] - posterior.draws.party.list[["JSP_HOR"]][, 2:38]), 2)

# average distance between the LDP and the DSP (HOR)
round(mean(posterior.draws.party.list[["LDP_HOR"]][, 2:38] - posterior.draws.party.list[["DSP_HOR"]]), 2)

# probability that the JCP is more left in 1959 than in 1989 (HOC)
round(mean(posterior.draws.party.list[["JCP_HOC"]][, "1959"] > posterior.draws.party.list[["JCP_HOC"]][, "1989"]), 3)

# probability that Komeito is more right in 2019 than in 2012 (HOR)
round(mean(posterior.draws.party.list[["Komeito_HOR"]][, "2019"] > posterior.draws.party.list[["Komeito_HOR"]][, "2012"]), 3)

party.estimates <- c()
for (i in 1:length(party.names)) {
  HOR.id <- which(party.levels == paste0(party.names[i], "_HOR"))
  HOC.id <- which(party.levels == paste0(party.names[i], "_HOC"))
  if (length(HOR.id) == 1) {
    party.estimates <- rbind(party.estimates, 
                             data.frame(party = paste0(party.names[i], "_HOR"), 
                                        year = c(party.begin[HOR.id]:party.end[HOR.id]) + 1958, 
                                        point = point.estimates.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOR.id, "\\,"))], 
                                        lower = credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOR.id, "\\,")), 1], 
                                        upper = credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOR.id, "\\,")), 2]))
  }
  if (length(HOC.id) == 1) {
    party.estimates <- rbind(party.estimates, 
                             data.frame(party = paste0(party.names[i], "_HOC"), 
                                        year = c(party.begin[HOC.id]:party.end[HOC.id]) + 1958, 
                                        point = point.estimates.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOC.id, "\\,"))], 
                                        lower = credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOC.id, "\\,")), 1], 
                                        upper = credible.intervals.party[stri_detect_regex(names(point.estimates.party), paste0("theta\\[", HOC.id, "\\,")), 2]))
  }
}
write.csv(party.estimates, file = "3-a_party_estimates.csv", row.names = FALSE)

#### committee loadings ####
posterior.draws.committee <- postprocessing.out[, stri_detect_regex(colnames(postprocessing.out), "beta")]

# posterior means and 95% HPD intervals for committee loadings
point.estimates.committee <- colMeans(posterior.draws.committee)
credible.intervals.committee <- matrix(NA, ncol(posterior.draws.committee), 2)
rownames(credible.intervals.committee) <- colnames(posterior.draws.committee)
for (i in 1:ncol(posterior.draws.committee)) {
  credible.intervals.committee[i, ] <- HPDinterval(mcmc(posterior.draws.committee[, i]))
}

# Figure 3: committee loadings over time (means and 95% HPD intervals)
committee.labels <- c("Cabinet", "General Affairs", "Judicial Affairs", 
                      "Foreign Affairs and Defence", 
                      "Financial Affairs", "Education, Culture, and Science", 
                      "Health, Welfare, and Labour", 
                      "Agriculture, Forestry, and Fisheries", 
                      "Economy and Industry", "Land and Transport", 
                      "Trans, Info, Land, and Envir", "Environment")

png("Figure_3.png", 
    width = 6, height = 4.57, unit = "in", pointsize = 8, res = 2000)
layout(matrix(1:12, 3, 4, byrow = TRUE))
par(mar = c(1.5, 1.5, 2, 0.2), lwd = 0.5)
for (i in 1:length(committee.labels)) {
  plot(NULL, NULL, type = "n", xlim = c(1, 61), ylim = c(0, 0.75), 
       main = committee.labels[i], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = seq(0, 0.75, 0.25), lty = 3, col = "gray")
  polygon(c(committee.begin[i]:committee.end[i], committee.end[i]:committee.begin[i]), 
          c(credible.intervals.committee[stri_detect_regex(names(point.estimates.committee), paste0("beta\\[", i, "\\,")), 1], 
            rev(credible.intervals.committee[stri_detect_regex(names(point.estimates.committee), paste0("beta\\[", i, "\\,")), 2])), 
          border = NA, col = "#80808030")
  lines(committee.begin[i]:committee.end[i], 
        point.estimates.committee[stri_detect_regex(names(point.estimates.committee), paste0("beta\\[", i, "\\,"))])
  axis(1, at = seq(1, 61, 10), labels = NA, lwd = 0.5)
  mtext(c(seq(59, 99, 10), "09", "19"), side = 1, line = 0.5, at = seq(1, 61, 10), cex = 0.6)
  axis(2, at = seq(0, 0.75, 0.25), labels = NA, lwd = 0.5)
  mtext(seq(0, 0.75, 0.25), side = 2, line = 0.5, at = seq(0, 0.75, 0.25), cex = 0.6)
}
dev.off()

committee.short <- c("cabinet", "general", "judicial", "foreign", "financial", 
                     "education", "health", "agriculture", "economy", "land", 
                     "land.2", "environment")
posterior.draws.committee.list <- list()
for (i in 1:length(committee.short)) {
  posterior.draws.committee.list[[i]] <- posterior.draws.committee[, stri_detect_regex(colnames(posterior.draws.committee), paste0("beta\\[", i, "\\,"))]
  rownames(posterior.draws.committee.list[[i]]) <- NULL
  colnames(posterior.draws.committee.list[[i]]) <- seq(committee.begin[i], committee.end[i]) + 1958
  names(posterior.draws.committee.list)[[i]] <- committee.short[i]
}

# average loading of the Economy and Industry Committee
round(mean(posterior.draws.committee.list[["economy"]]), 2)

# average loading of the Land and Transport Committee
round(mean(posterior.draws.committee.list[["land"]]), 2)

# average loading of the Foreign Affairs and Defence Committee
round(mean(posterior.draws.committee.list[["foreign"]]), 2)

# probability that the Cabinet Committee's loading is higher in 2019 than in 2000
round(mean(posterior.draws.committee.list[["cabinet"]][, "2000"] < 
             posterior.draws.committee.list[["cabinet"]][, "2019"]), 3)

# probability that the Education Committee's loading is higher in 2019 than in 1998
round(mean(posterior.draws.committee.list[["education"]][, "1998"] < 
             posterior.draws.committee.list[["education"]][, "2019"]), 3)

# committee ranking probabilities by year
committee.rank.results <- array(NA, c(61, 11, 3))
for (i in 1:61) {
  if (i < 22) {
    committee.draw.year.i <- posterior.draws.committee[, paste0("std.beta.beta[", 1:10, ",", i, "]")]
    committee.rank.year.i <- 11 - t(apply(committee.draw.year.i, 1, rank))
    committee.rank.results[i, 1:10, 1] <- apply(committee.rank.year.i, 2, function(x) mean(x == 1))
    committee.rank.results[i, 1:10, 2] <- apply(committee.rank.year.i, 2, function(x) mean(x < 4))
    committee.rank.results[i, 1:10, 3] <- apply(committee.rank.year.i, 2, function(x) mean(x < 7))
  } else {
    committee.draw.year.i <- posterior.draws.committee[, paste0("std.beta.beta[", c(1:10, 12), ",", i, "]")]
    committee.rank.year.i <- 12 - t(apply(committee.draw.year.i, 1, rank))
    committee.rank.results[i, , 1] <- apply(committee.rank.year.i, 2, function(x) mean(x == 1))
    committee.rank.results[i, , 2] <- apply(committee.rank.year.i, 2, function(x) mean(x < 4))
    committee.rank.results[i, , 3] <- apply(committee.rank.year.i, 2, function(x) mean(x < 7))
  }
}

# Figure A4: committee ranking probabilities by year (top-1, top-3, top-6)
committee.labels.2 <- c("Cabinet", "General Affairs", "Judicial Affairs", 
                        "Foreign Affairs\nand Defence", 
                        "Financial Affairs", "Education, Culture,\nand Science", 
                        "Health, Welfare,\nand Labour", 
                        "Agriculture, Forestry,\nand Fisheries", 
                        "Economy and Industry", "Land and Transport", 
                        "Environment")

png("Figure_A4.png", 
    width = 6, height = 3, unit = "in", pointsize = 8, res = 2000)
layout(matrix(1:12, 2, 6, byrow = TRUE))
par(mar = c(1.5, 1.5, 3, 0.5), lwd = 0.5)
for (i in 1:11) {
  plot(NULL, NULL, type = "n", xlim = c(1, 61), ylim = c(0, 1), 
       main = committee.labels.2[i], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = seq(0, 1, 0.2), lty = 3, col = "lightgray")
  abline(v = 22, lty = 3, col = "lightgray")
  lines(1:61, committee.rank.results[, i, 1])
  lines(1:61, committee.rank.results[, i, 2], lty = 2)
  lines(1:61, committee.rank.results[, i, 3], lty = 4)
  axis(1, at = seq(1, 61, 10), labels = NA, lwd = 0.5)
  mtext(c(seq(59, 99, 10), "09", "19"), side = 1, line = 0.5, at = seq(1, 61, 10), cex = 0.6)
  axis(2, at = seq(0, 1, 0.2), labels = NA, lwd = 0.5)
  mtext(seq(0, 1, 0.2), side = 2, line = 0.5, at = seq(0, 1, 0.2), cex = 0.6)
}
par(mar = c(0, 0, 0, 0), pty = "m")
plot.new()
legend("center", legend = c("1st", "3rd or higher", "6th or higher"),
       lty = c(1, 2, 4), bty = "n", cex = 1.2)
dev.off()

#### assessing uncertainty in LDP and JCP ideal-point shifts (1996-2019) ####
# validation-related check added here to assess whether apparent negative correlations
# for the LDP and JCP reflect genuine shifts or posterior uncertainty
LDP.draw.Kato.period <- posterior.draws.party.list[["LDP_HOR"]][, as.character(c(1996:2019))]
round(mean(LDP.draw.Kato.period[, which.min(colMeans(LDP.draw.Kato.period))] < 
             LDP.draw.Kato.period[, which.max(colMeans(LDP.draw.Kato.period))]), 3)

JCP.draw.Kato.period <- posterior.draws.party.list[["JCP_HOR"]][, as.character(c(1996:2019))]
round(mean(JCP.draw.Kato.period[, which.min(colMeans(JCP.draw.Kato.period))] < 
             JCP.draw.Kato.period[, which.max(colMeans(JCP.draw.Kato.period))]), 3)
